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A HYDRODYNAMIC LIMIT FOR CHEMOTAXIS IN A GIVEN 
HETEROGENEOUS ENVIRONMENT 

IN HONOR OF WILLI JAGER, ON THE OCCASION OF HIS 75TH BIRTHDAY. 


STEFAN GROSSKINSKY, DANIEL MARAHRENS, AND ANGELA STEVENS 


Abstract. In this paper the first equation within a class of well known chemo- 
taxis systems is derived as a hydrodynamic limit from a stochastic interacting 
many particle system on the lattice. The cells are assumed to interact with 
attractive chemical molecules on a finite number of lattice sites, but they only 
directly interact among themselves on the same lattice site. The chemical envi¬ 
ronment is assumed to be stationary with a slowly varying mean, which results 
in a non-trivial macroscopic chemotaxis equation for the cells. Methodologi¬ 
cally the limiting procedure and its proofs are based on results by Koukkus 
m and Kipnis/Landim d- Numerical simulations extend and illustrate the 
theoretical findings. 
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1. Introduction and result 

In this paper we derive chemotaxis-like equations as a hydrodynamic limit of a 
stochastic lattice gas. Chemotaxis describes the directed motion of mainly biolog¬ 
ical species towards higher or lower concentrations of chemical signals. Here we 
consider positive chemotaxis of cells , i.e. motion towards higher concentrations of 
a chemical signal, therefore the chemical signal is denoted as chemo-attract ant 

Keller and Segel proposed in [16j a phenomenological chemotaxis model on the 
macroscopic level for the aggregation and self-organization of the cellular slime 
mold amoeba Dictyostelium discoideum (Dd). This can be written as 

(1.1) d t p = V- (fcGo^Vp-xO^W) 

= v • ( k(p , i?y vp - x(p, , 

(1.2) d t d = Ad + r(p, $), 


Date : August 17, 2015. 

Key words and phrases, chemotaxis, interacting stochastic many particle system, hydrody¬ 
namic limit, stochastic lattice gas, block estimates. 


1 



2 


S. GROSSKINSKY, D. MARAHRENS, AND A. STEVENS 


where p is the (volume) density of amoebae, P is the density of the attractive chem¬ 
ical molecules, k and X = x/p are functional parameters describing the strength 
of random motion of the cells and their chemotactic sensitivity, respectively, and 
r comprises the reaction mechanisms. In [16], this cross-diffusion system was mo¬ 
tivated by the macroscopic phenomenon being experimentally observed, namely 
movement of the amoebae in direction of higher concentrations of a chemical sig¬ 
nal, i.e. their movement up chemical gradients. Such general types of chemotaxis 
systems are relevant also in the context of other biological species. One expects 
this system to be an accurate description for chemotaxis phenomena occuring in 
systems of many cells and signal molecules. There is a large literature on the formal 
derivation of macroscopic equations from microscopic particle systems in such set¬ 
tings. Rigorously this has been achieved in much fewer cases, e.g. by hydrodynamic 
or so-called moderate limits. The lack of an ellipticity condition for the limiting 
PDE-system, respectively the strong clustering of the cells is a major technical 
problem in this case. 

The first rigorous derivation of a general class of chemotaxis systems from a sto¬ 
chastic interacting many particle system was given in £3, where the cells and the 
chemical molecules interact moderately with each other in the sense described in 
[28] . The motion of the cells and the molecules are governed by interacting sto¬ 
chastic differential equations, and production and decay of the chemical molecules 
are modeled via Poisson point processes. In explicit error estimates could be 
given, which show that the total number of particles does not have to be too large 
for the limiting PDE-system to be a good approximation for the particle model. The 
moderate rescaling requires that in the limiting procedure, when the total particle 
number (cells and chemical molecules) tends to infinity, the main range of inter¬ 
action of the particles tends to zero. The number of particles sensed in this main 
range of interaction tends to infinity too, while this number is still only a vanishing 
fraction of the total particle number. This is sufficient for correlations to become 
small enough in the limit. The main technical complication in the proof is the 
cross-diffusion structure of the chemotaxis system. The ellipticiy condition, heavily 
used in [28], is no longer valid in this case. To overcome this problem, two shadow 
systems were introdcued, in order to freeze the critical non-linearity and a priori es¬ 
timates were derived. This technique applies also to more general systems of PDEs. 


In principle it would be desirable to also derive chemotaxis like systems as hydro- 
dynamic limits in the sense of m- In this paper, we rigorously derive equation 
(1.1) from a microscopic stochastic many particle system on a lattice, so only one 
equation, not the full system. The cells are assumed to interact via a finite number 
of lattice sites with the chemical molecules. This corresponds to the assumption 
that the cells detect chemical molecules in a fixed small region around themselves. 
They do not change the chemical environment in this case. The interaction of the 
cells among themselves takes place only at the same lattice site and not among 
cells located on neighboring sites. This interaction is described by a function g 
which among others fulfills conditions comparable to an ellipticity condition for 
the related limiting PDE. These will be specified later. So we consider the case, 
where a too strong clustering/aggregation of cells - a typical and important effect 
of chemotaxis and self-organization in Dd - is avoided. 

Such an interaction with a finite number of lattice sites introduces specific mathe¬ 
matical difficulties. One has to deduce a weak equation for the limit density from 
empirical measures (see Section [2| which converge weakly, cf. the discussion in 
[27, Introduction]. To be more precise, we are modeling the cells/particles by an 
interacting particle system as introduced in [34] and derive limit equations via a 
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hydrodynamic limit, analogous to the procedure in PH- Interacting particle sys¬ 
tems in this sense are continuous-time Markov jump processes which also involve 
discrete particles moving on a lattice. 


So far we were not able to derive the full system (1.1), (1.2) as hydrodynamic limit. 
It seems to be a major challenge to identify suitable particle models with invariant 
product measures. This would be one option in order to prove a hydrodynamic 
limit for systems of equations such as 0 and ( |1.2| ). Hydrodynamic limits without 
product measures for single-species monotone particle systems where established in 
H 0. Stochastically monotone systems preserve a partial order on state space 
over time, which allows the use of coupling techniques following first results in [32] . 
Monotonicity usually leads to homogeneous mass distributions, and indeed it has 
been shown recently that particle systems with invariant product measures that 
exhibit a particular form of clustering are necessarily non-monotone m • 

A generic example which exhibits invariant product measures is given by the zero- 
range process (ZRP), where the jump rate of particles at any given site depends 
only on the occupation number at this site. It is known, see [13], that multi-species 
ZRPs have invariant product measures if the jump rates of particles satisfy certain 
symmetry relations, i.e. the Onsager relations. However, in our case this would 
not allow for purely diffusive motion of the chemo-attract ant in 0. Instead, the 
attractive chemical molecules would need to undergo a kind of chemotactic motion 
too, in order for the Onsager relations to hold. From the modeling point of view 
this is definitely not the case here. 

Another method to derive a system of a PDE coupled to an ODE via a hydrody¬ 
namic limit without the use of invariant product measures was introduced in m 
for a specific two-species cellular system on the lattice. The authors new method 
employs energy estimates for the mesoscopic empirical averages of the particles oc¬ 
cupation number by which they obtain H\ a-priori bounds and thus are able to 
derive a substitute for the two block estimates. For the most crucial terms in their 
setting homogenization techniques are used, which play the role of the one block 
estimates. For a discussion of the respective particle model itself see also [23]. It 
would be interesting to see whether in our case one can find scaling regimes where 
the existence of gradients (and not only densities as in [TO]) can rigorously be proved 
for the full limiting PDE system. 

One typical approximation of system 0,0 is to assume a quasisteady chemo¬ 
attractant density, i.e. to set the left hand side of 0 equal to zero. Formally, this 
is justified by assuming that the diffusivity of the chemo-attract ant is much larger 
than the random motion of the cells, see e.g. [14 . There are also recent results on 
random walks in dynamic random environments driven by exclusion processes [3], 
which are based on separation of time-scales argument. Otherwise there are only 
very few rigorous results on hydrodynamic limits for two-species systems, see e.g. 
[ 117 ], since such processes are in general also non-monotone. 

It is obvious from a modeling point of view, how a full chemotaxis model for the 
behavior of particles on a lattice could look like, compare e.g. [29 . One ansatz 
relates to attractive reinforced random walks for many particles with diffusion and 
decay of the attractive weight. Deriving a hydrodynamic limit for such, or similar 
models would be the final goal. We expect the Keller-Segel model to also hold in 
this context, but a proof is still missing. 


As a first step towards understanding chemotaxis as a hydrodynamic limit, we 
therefore consider chemotaxis of cells in a stationary, but random, environment. In 
order to obtain a non-trivial macroscopic chemotactic motion, we impose a slowly 
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varying mean on the random chemical environment. There have been several stud¬ 
ies of hydrodynamic limits of particles in random media: The hydrodynamic limit 
for a ZRP a with stationary, ergodic environment on the sites has been obtained in 
[18], and corresponding large deviations have been considered in [19] . The case of 
a stationary, ergodic environment on the edges was treated in m- The hydrody¬ 
namic limit for the ZRP with slowly varying but not random environment is given 
in |9j. For a result unifying all these situations of “locally convergent” media in the 
special case of the simple exclusion processes compare [15] . Our main result in the 
present paper concerns the modeling of chemotaxis via a ZRP in a random envi¬ 
ronment with slowly varying mean and it is summarized in Theorem |2.1[ Section 
HI A special case of this result is given in Corollary [hi] below. 


At each site x G on the periodic lattice = {1 ,...,N} d one distributes 
((x) G N molecules of the chemo-attractor such that ((x ) is a Poisson random 
variable with parameter 'd(^) and these are kept fixed for all time. Let / : N -A [a, b] 
for some a, b > 0 and set 


/(<?(#)) :=E 


/(C(*))J 


-1 


Now we specify the type of chemotactic motion we consider here for the cells. We 
start with an initial distribution of rj(x) cells at x G T^. All cells perform inde¬ 
pendent random walks with a site-dependent jump-rate which is N 2 f(((x)). Thus 
each cell remains at its current site x G for an exponentially distributed random 
waiting time with parameter N 2 f(((x)) and then jumps to a random neighboring 
site on the lattice. In this way, particles are more likely to stay at a site x if f(((x)) 
is large and therefore / is one way of formulating the microscopic behavior which 
results in a chemotactic drift towards higher concentrations of the chemical signal. 
To take the limit as N -A oo, we embed the discrete torus into the continuous 
torus T d = M. d /Z d via r G Thus the number of cells in an interval I C T d is 
given by 

#{cells in 1} = g(x), 

xeNinT% 

where NI is / rescaled by N. The cell density in I is then obtained as 
E x eNiV(x)/(\I\N d ). 


Corollary 1.1. Let r\ q(x) be distributed with a profile po(u), i.e. for each interval 
I C T d , the density converges in probability: 


lim 

AT—>• oo 


#{cells in 1} 
|/|7V d 


Sr ( 


u) du. 


Then for all later times t > 0 and all intervals I, almost surely with respect to the 
distribution of chemical molecules it holds that 


lim 

AT—>> oo 


#{cells in 1} 
\I\N d 


= JjPtiu) 


du 


in probability, where p t solves 

(1.3) d tPt (u) = A[f(#(u))p t (u)] = V • [/(t?(u))V/Jt(«) + /'(0(u))pt(u)V0(u) 
which is of the form <o>- 


For the derivation of this result from Theorem 12. II we refer to the discussion in 
Section [4] Let us first make several remarks concerning this result. 
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Remark 1.2. (1) The assumption of independent random walks treated in Corol¬ 
lary [TTl] yields a linear limit equation and is much weaker than necessary, see The¬ 
orem |2.1| It is presented here in order to give the reader an intuition for the 
connection between the macroscopic bias /($( u )) induced by the chemical signal 
and the related microscopic bias /(£(#)). 

(2) Note that instead of the Poisson distribution, we can let the environment be 
given by any distribution in N depending continuously (in distributional topology) 
on its parameter. 

(3) In case of a stationary chemical environment, the solution is always global and 
no blow-up can occur, thus excluding a prominent feature of the full chemotaxis 
system, where the possibility of finite time blow-up in two space dimensions for a 
suitable parameter range is of important biological relevance for self-organization 
phenomena in the cellular slime mold amoebae Dd, i, mi. 

(4) Finally, let us mention that some microscopic descriptions for chemotactic cell 
motion are based on the feature that the cells can detect concentration gradients 
along their cell surface. Our result supports the idea that chemotactic effects can 
also occur if the cells only sense the absolute values of the concentration of the 
chemo-attractant without needing to sense its (local) gradients. An alternative ap¬ 
proach would be to let each particle perform asymmetric random walks (and its 
nonlinear versions involving a zero range interaction #(•), see below) with jump rates 
across edges of the lattice proportional to p(x,y). Here x, y are neighboring sites, 
and p(x,y) is a random variable, which in our case is determined by the attractive 
chemical environment £(#)., see also [29] and the microscopic approaches in [2]. A 
typical example for an asymmetric random walk would be p(x, y) = f(((y) — ((x )), 
representing a microscopic gradient for an appropriate function /. This approach 
would lead us to consider problems of the type investigated in [12.. 

(5) For positive chemotaxis effects, i.e. clustering/aggregation of the cells the sign 

of the second term on the right hand side of © is crucial. Further, already 

in [33] h was pointed out that the diffusivity /, k and the chemotactic sensitivity 
—/', x °f the cells are not comp letely independent functionals. A suggestion in [33] 
is that the diffusivity / in ( |l.3| relates to the chemotactic sensitivity — /', which is 
supposed to have a positive sign in case of positive chemotaxis, as follows: 

/' = —xo f&($) for a suitable function 4>, e.g. $($) = $ and xo > 0. 

This results in / = Cexp(—xo^C#)), e.g. / = C exp(—Xo'#)- 
So for large chemical concentrations / becomes small. 

The outline of the paper is as follows. In Section [2] we present a general result 
on hydrodynamic limits in a random environment with slowly varying mean. Its 
proof is given in Section [3] where we highlight the differences to the proofs found in 
[Mill]. Then in Section^) we show how to deduce Corollary © from the general 
result and present numerical simulations of the particle system. 


2. The general result 


In this paper, we apply the entropy method as given in [25] to obtain a generalization 
of equation (1.3), namely for the chemotactically moving cells, in a stationary, yet 


random, environment, which represents the density of a heterogeneously distributed 
attractive chemical signal. To avoid technicalities concerning boundary conditions, 
we shall work in a periodic domain. Thus we consider a Markov (Feller) process, 
describing the motion of the cells, with state space N TiV , where Tjy = {1,..., N} d 
denotes the perodic lattice with N + 1 = 1. Elements of the state space will 
be called particle configurations and denoted by the Greek letters 77 , £, £. Thus 
rj(x) G N denotes the number of particles of the configuration 77 at site x. The 
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elements of will be denoted by the letters x,y,z and are called microscopic 
variables. Two microscopic sites x, y G are called neighbors , in short x ~ y, if 
\x — y\ = 1. The discrete torus is embedded in the continuous torus T d = R d /Z d 
via x i— y x/N. Elements of the continuous torus will be denoted by u and are called 
macroscopic variables. For simplicity, we consider only symmetric nearest-neighbor 
jumps; then the Markov process is given by the generator 

(2- 1 ) £jv/0) = E 9^ x ))Px (f(v X ’ V ) - f(v)) 

x~y 

for all / G Cfr(N T ^), where the sum is taken over all (ordered) pairs of neighbors 
x and y. Here y x,y denotes the configuration obtained from rj after one particle 
has jumped from site x to y and p£ > 0 describes the chemical environment. 
Throughout this article, we assume that 

P% -v(§) + q x > 

where v G C 1 (T d ,M) describes the slowly varying mean of the environment and 
(q x )xez d C is a uniformly bounded stationary and ergodic sequence of random 
variables with zero mean. In order to avoid degeneracies, we assume strictly positive 
jump rates, i.e. without loss of generality we assume p™, q x E [a, b] with 0 < a < b. 
We denote the law of q by m, so that m is a probability measure on R z . 

It is also possible to think of p£ as a random variable. Note, however, that for 
technical reasons in what follows we fix the ergodic part q of the environment 
independently of N. The dependence of the jump rates on the number of cells is 
given by the function g : N -A [0, oo). In addition to the standard condition 

g{n) = 0 n = 0 

to avoid degeneracies, we make the following regularity assumptions on g through¬ 
out the paper. 


Assumption 1. (i) Suppose that g is uniformly Lipschitz-continuous, i.e. there 
exists a constant g* such that 


sup|sr(n + l) - g{n)\ < g*. 

n£ N 


(ii) Further assume that g grows at least linearly, i.e. there exists go > 0 such that 


inf 

n£ N 


gjn) 

n 


> go- 


These assumptions are not optimal but they are standard in the literature, see e.g. 
m Theorem 5.1.1], on which our proof is based. Under these assumptions, the 
process with generator <H3 has invariant product measures v N,p that satisfy the 
detailed balance condition 


Px 9(n)v N ’ p [g{x) = n\u N ’ p [r](y ) = k] 

= Py g(k + 1 )v N ' p [g{x) = n - l]u N,p [i](y) = k + 1] 
for all neighbors x, y G and k,n G N. 

Indeed for any (p > 0, there exists such a measure given by 
<2 ' 2) 

x£T% 

where g(n)\ = g(l)g(2)... g(n ), g(0)! = 1 and 


»)~M V{X) 

)~W)g(v (*))! 


(2.3) 


z(<p) = E 




q(n)\ 

n=0 v J 
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is the partition function. The detailed balance condition ensures that the generator 
C N in (2.1) is L 2 {y^ p )~ selfadjoint and in particular the measure n^ ,p is an invariant 
product measure. Let zA denote the one-site marginal without environment, i.e. 


d( n ) = 




Z(<p)g(n)\ 


The associated density then is 

(2.4) 


M(<p) mEM 0)]. 


Assumption [I] implies that Z is finite on [0, oo). Hence the product measure v^ p 
exists for all cp G [0, oo) and environments (Px) x e t^- The parameter ip > 0 is called 
fugacity, and it controls the expected particle density which is given by 


(2.5) 


R(u,<p) :=E m [M(^)^ 


for all u G T d . The family of distributions (2.2) is also called the grand-canonical 
ensemble. One can show that for each u G T d , the function R(u, •) is strictly 
increasing, cf. m, and we then define <f>(w, •) to be its inverse function. Note that 


E„ 


K?* ,D?(°)] 


p and E m Ew, P [g(ri(x))p*] =$(§,<p). 


We shall see that the limit equation is given by 

(2.6) d t p t {u) = A ®(u,p t (u)). 


Finally we need a suitable notion of distance between probability measures, which 
in our case is given by the relative entropy. Consider measures u, p G P(N Tn ) such 
that (i is absolutely continuous with respect to v (i.e. p <C v). Then the relative 
entropy is given by 

H(g\v)= f log^(r?) dp(ri), 

JWn nU 

where ^ denotes the Radon-Nikodym derivative of p with respect to v. The main 
result of this section is the following. 


Theorem 2.1. Let p{f G P(N t ^) be the initial datum of the particle process with 
an associated initial profile po G L°°(T d ), i.e. it holds 


lim Pq 

N^oo 


fd 52 G (§)r](x) - f G(u)p 0 (u) du 
xeTt Jjd 


>6 =0 


m-almost surely 


for all G G C(T d ) and S > 0. Furthermore we suppose the bounds 


H(p*\v^)<CN d 


and 




52 d{x) 2 

l xer d N 


< CN d 


to hold for some (and hence all) > 0. Denote by G P(N T ^) the measure 
obtained from the evolution of the ZRP with rate function g and by pt € L°°(T d ) 
the density obtained from equation (|2.6| . Then under Assumption^ it holds that 


Jim T‘t 

N—> OO 


N 


1 

w 


E 


- [ G(u)p t (u) du 

J Jd 


> 6)=0 


m-almost surely 


for all t > 0, G G C(T d ) and S > 0. 
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Remark 2.2. In order to fully connect Corollary |1.1| with Theorem 2.1, we note 
that weak convergence of the empirical measures implies convergence of the densities 
over intervals by the Portmanteau theorem [6j Theorem 2.1]. This follows from the 
absolute continuity (with respect to the Lebesgue measure) of the limit measure 
pt(u) du. 


To prove this result, we adapt the method in m to a random environment with 
slowly varying mean. The paper RBI is based on the entropy method given in [25 
and proves a hydrodynamic limit for a zero range process in a stationary (ergodic) 
random environment. The case of a zero range process in a slowly varying deter¬ 
ministic environment was treated in [9]. Thus, from a purely technical point of 
view, our result is a combination of the previous two. Indeed, our method is close 
to the methods presented in m • Therefore we mainly highlight the differences 
here. The basic idea of the proof lies in the fact that the slowly varying mean is 
locally almost constant and hence in small boxes we are essentially in the situation 
of PS]. 


3. Proof of Theorem 1271 


The first step of the proof is to obtain a priori bounds and tightness of the particle 
process. Indeed we can understand the convergence result of Theorem |2.1| in terms 
of the empirical measure 


a v ( du ) = jfd 57 e M +’ 


ceTl 


where S u denotes a Dirac delta, u E T d , and M.+ is the space of all positive Radon 
measures on T d . For allt > 0, the configuration ry is a rand om v ariable distributed 
according to . Note that it suffices to prove Theorem 2.1 for all 0 < t < T 
with a fixed, but arbitrary, T > 0. Since for any environment p the ZRP is a 
jump process on the state space, (rjt)te[o,T] is a random variable on the path space 
D([0,T];N Td ). Here D([0,T];N Td ) denotes the space of functions [0,T] -A N jd 
that are right-continuous in time with left limits. We denote the distribution of 
{Vt)te[o,T] by p N . In the same vein, (&% t )te[o,T] is a random variable on the path 
space D([0, T]; jM+); let Q N,P denote its distribution. The entropy dissipation of 
/if can be defined in terms of the Dirichlet form 




? - 


): L 


,N 


d ^ dv N *. 


dv L 


N,p 


Throughout this paper we suppose that the assumptions of Theorem 2.1 are satis¬ 
fied. 


Lemma 3.1. It holds that 

H(fi? |zA’ p ) < CN d for all t € [0, T] and I [ V(^\v^’ p ) dt < CN d . 

T J o 

Lemma 3.2. For any environment p, the sequence of probability measures ( Q N,P ) tvgn 
is tight. 

Lemma 3.3. For any environment p, all limit points Q* ,p of (Q N,p )NeN are sup¬ 
ported on {7r G D([0, T]; *Ad+(T d )) : 7 r t (du) < du}, i.e. at each time t, the limit 
measures are absolutely-continuous with respect to the Lebesgue measure. 

Lemma 3.4. For any environment p, if 7r t (du) = pt(u ) du is distributed according 
to Q* ,P ; then p t (u ) is the unique solution of ( |2.6| ) in L 2 ((0,T) x T d ). 
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The proofs of Lemma |3.1| - |3.4| are similar to the proofs of the corresponding 
results for the usual ZRP, see Lemma 2.1 - 2.4 in [18] and Lemma V.1.5 and V.1.6 


in na. Once these results have been obtained, Theorem |2.1| is a straight-forward 
consequence. The main difficulty consists in the proof of the so-called replacement 


lemma, which is the key to Lemma 3.4 In order to state the replacement lemma, 
we need some additional notation to describe averages over mesoscopic blocks. Set 

1 V 


and 


Lemma 3.5 

that 


n{x) 


v xM-= jytLw E Py9(v(y))-$(§,v l (x)) 


(2l + l) d ^ 


y{y) 


\x-y\<l 

(Replacement Lemma). For every S > 0, 


m-almost surely it holds 


lim sup lim sup /i 

e—^0 TV—>- oo 


N 


(jf 


l 

W 


E 


V x,cnM dt> 6) = 0, 


where t x denotes the translation by x £ T(( 


N ’ 


In the remainder of this section we will mainly prove Lemma 3.5, where we sketch 
the modifications necessary to adapt the proof in m in order to take into account 
the slowly varying average of the environment First we prove the replacement 

over small boxes of size l G N. 


3.1. One block estimate. 


Lemma 3.6 (One Block Estimate). It holds that 


lim sup lim sup E^ 

Z—»• oo AT—>• oo 


I w E v xM) dt 

xefl J 


= 0 


m-almost surely. 
Setting 


fl(v) = ywjiv), pN (v) = f ft N (v ) dt, 

,jP Jo 


we see that 


xeTt xeT<* 


E fJ 


Lemma [3.1[ together with convexity of the entropy and Dirichlet form, yields the 
bounds 

H n ’P ( F n ) := H ( F n djy N , P | du N, P ^ < CN d ? 

V N 'P(F n ) : = V(F n dis*’ p \ dis*’ p ) < CN d ~ 2 . 

For future reference, we also note the following expression of the Dirichlet form via 
the density F N . It holds that 


(3-1) 


(3.2) x> n,p (F n ) = Y, L. \Px9(jl{x))[J f»(tT’V) - JfFI)] 2 dv F . 

X ~y J ^ N Z 


The next lemma allows us to restrict ourselves to bounded particle configurations. 
Its proof is analogous to the proof of m Lemma V.4.2] with the entropy inequality 
as additional ingredient. 
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Lemma 3.7. Under the conditions of Theorem \2. 1\ we have 
lim sup lim sup lim sup sup / F V V^i(r])x{r,‘(x)>A}(v)f(v) dA’ p (r]) = 0- 

A—>-00 l —yoo AT—>■ 00 p,f Jw^n d 

X(z Ji jy 

where the supremum is taken over all densities f such that H N,P (F N ) < CN d , 
X) n >p(F n ) < CN d ~ 2 and all environments p = C [a, b] TN . For every set A, we 
denote by xa the characteristic function of A. 


Therefore we just need to prove 

lim sup lim sup -T ^ f V p lA (v)F N (v) dv*’ p {ri) = 0 

l — v 00 V—>-oo d Jn iv 

ccGN T -ZV 

m- almost surely, where we set Vf t A (rf) := V^ l (rf)X{ r fi(x)<A}{ r i)- I n fact, due to the 
bounds (3.1) on the Dirichlet form of F N , it suffices to find an upper bound for 

(3-3) W E J N< V *i,aW fN (v) dv^(v) - 7 CN*- d V N *(F N ) 


for all 7 > 0. Next we restrict the problem to translations of the small box A i := 
{—21 ,..., 21}. Denote by the A x y -marginal of u^ p and by F x j(rj) the density 

of the -marginal of the measure F N (rj)du^ ,p ( 77 ) with respect to v* ,l,p . We will 
sometimes drop the index for x = 0 for such quantities when taken at the origin. 
Furthermore we define a Dirichlet form on A x j by 


£*£,*(/) = Y tyM)’ where we have set 

y~Z 

y,zeA Xj i 

dyM) = \ jf Py giviv)) (Vf(v v,x ) - v7b)) dv^ p . 

By convexity of the ”bond”-Dirichlet forms I p y , we have I p , z (F x j) < I P , Z (F N ) for 
all neighbors x,y G A x j and thus 

xeT% 

Hence instead of expression ( |3.3| |, we just need to find an upper bound for 

( 3 - 4 ) W E { d K’ l ’ p W - 7 

Next we need to take care of the random part of the environment. The method 
to keep track of the environment is due to [18] and constitutes the main deviation 
from the proof of the hydrodynamic limit for the usual ZRP as exhibited in m- 
Thus we fix a, S > 0 and let n G N be sufficiently large such that ^ < S. Divide the 
interval [a, b\ into sub-intervals of length not greater than S(b—a ) via I d = [f3j, /3j+i) 
for 0 < j < n — 2 where 

/3j = a + (b- a)^ (j = 0 , 1) 

and I^i = [/? n _i, &]. Fix k < l and let L = [(21 + 1)/(2 k + l )] d , where [x\ denotes 
the Gaussian bracket, i.e. the largest integer smaller than or equal to x G M. Now we 
divide A/ into disjoint cubes of the form x + A&, where we take Bi, i = 1,..., L — 1, 
such that 

Bi C A i, Bi n Bj = $ for i ^ j , and Bi = Xi + A/ for some 27 G Z d . 
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Then set Bl = A i \ [jf =1 1 Bi and without restriction take x\ = 0, i.e. B\ = A&. 
Also set Bi(x) = x B{ for x G Z d . For x G T^-, a G (0,1) and q G [a, 6] Az , set 

1 


:= fOL 1 53 


(2fc + l) d 


zEBi(x) 


(3.5) 


eA, : | N l x ’j,t(q) ~ Em[X/«(g*)] < a, j = 0,... ,n - l|, 

A x k a '■= j(g*)*€Ai : J^J2 x A l ’ k M q ) > 1 - a| • 


Since Vf' A < C(A) is bounded, the expression (3.4) is bounded form above by 

(3.6) 

w E icm 2 KAF x ,i) 

+ COX 1 - x^Xg)) j- 

Now we take care of the non-random part of the enviroment. To this end, set 


Xl,k,8 


{ (Pz)ze A, : p z = Qz +v z ,q e A l ^ s , sup \v x+z - v(u )| < a} 

L zeAi J 


Since v is smooth, in particular (uniformly) continuous, for every a we can choose 
N large enough such that v x + z := v({x + z)/N) does not differ from v(x/N) by 


more than a. Thus (3.6) is bounded from above by 

(3.7) 


A E Ua^ (p)( L d^ l ’ p (v) -lCm 2 Ki( F ^ 

N xex- l KJNN 


+ C(A)(1 - X A ^(q))'j- 

Since the random part q of the environment is stationary and ergodic by assumption, 
it holds 

1 E i 1 ~ Xa«*'(9)) V rn(q i A 1 ^), 


N d 




which in turn vanishes by ergodicity in the limit as / -A oc and k 


oo. The 


remaining term in (3.7) is bounded from above by 


sup Xjif.*.* (p)( [ Jd V[uNU )A ^ F ^N]M dv^’ l ’ p (ri) 

ueT d [uN],u, a \J N T N L J> 

-icm 2 ^ uN] /F [uN]tl ) 

Again r_^ uN ]P is in A l Q k ^ 6 a for large enough N and the term involving <h(^Al,p) 
which appears in V^ uN ^ t A converges uniformly in u and p < A to <F(r, p). Thus we 
see that, up to an error that vanishes as N -A oo, it suffices to estimate 

(3.8) sup sup sup ([ VaO)/^) dv l f(rj) -^C{l)N 2 V p 0 ,(/)), 

uer d v eA l ’ k ’ s f wn A i J 

0 ,u,cx 
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where the inner supremum is taken over all densities / on N Az . For fixed l G N, the 
term Vq L a is only non-zero on the compact space of configurations in N A * of at most 
(21 + l) d A particles, the Dirichlet form is lower-semicontinuous and this property 
is conserved by the supremum. Hence the limit superior of (3.8) is bounded from 
above by 


(3.9) 


sup sup sup [ V t p A (v)f(v) dv l /(ri), 

ueT d vG A l ’ k ’ 5 f 

0,u,a. 


where now the inner supremum is taken over all densities / on N Az with vanishing 
Dirichlet form V p (/) = 0. A density / has vanishing Dirichlet form if and only if it 
is constant along all hyperplanes of a given number of particles. With this in mind, 
we introduce the canonical measures 


vfM ■) : = V ,P ( • I E*eA, v(x) = K) . 

Since the probability density / is constant on {^ xeAl rj(x) = K}, we can esti¬ 
mate (3.9) from above by 


(3.10) sup sup max 

u£peA l 0 ^ s a 0<AX(2Z+l) d A 


{Lit! 

- $ G >(2 wiy) dv iAv)\, 


where we have inserted the explicit form of Vq L a and applied the triangle inequality. 
For given u and p G we define p p K by 

U 11 ) \i ln ' ix)i = why- 

Note that 0 < p p K < C(A) is uniformly bounded in 0 < K < (21 + 1 ) d A and l G N. 
The next lemma concerns the closeness of the grand-canonical and the canonical 
measures. 


Lemma 3.8 (Equivalence of ensembles). For all F : N Afe -A M with finite second 
moments with respect to for all p < C(A), it holds that 


lim sup max 

l ^°°pe[a,b] zd °< K <( 2l + 1 ) dA 


K^[F(v)}-K r jFm 


0 . 


Proof. Thanks to the uniform bounds a < p x < 6, the equivalence of ensembles can 
be proved as explained in pTl Appendix 2], see also [I8j Lemma 3.3]. □ 

By applying Lemma |T8| to 


F(v) = 


1 


(2k + l) d 


E Px9(v(x)), 


x£Bi 


we just need to find an upper bound for 


(3.12) sup sup max 

u£T d N peA lA / a 0<^<(2Z+l) d A 


Lit 


(2k + l) d 

-$(l 


E Px 9 (v(x)) 


x£Bi 

K > 
(2l+l) d ) 


dv%(v) 
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The next two lemmas yield convergence of this expression in the limit as / -A oo 
and then k oo. In a certain sense, the first lemma describes the “ergodic” error 
between spatial averages and the quenched grand-canonical average. 

Lemma 3.9. For all cq S > 0, it holds that 


lim sup lim sup sup sup max 

k —>oo l —yoo uETfj p£A l Q ’ k,s 


{ L \ t / ipVip y ™ Mx)) - *i =°- 

The second lemma describes the difference between the quenched and annealed 
grand-canonical averages. 

Lemma 3.10. It holds that 


lim sup lim sup lim sup lim sup sup sup max 

<5—>-0 a—>-0 k —>-00 l—too pEAq k,S 0<-^E(2£+l) d A 




= o. 


Proof of Lemma \3~^ Under the law is l fp , the random variables (p x g(p(x))—p^ <: ) xe B i 
are independent with zero mean. Furthermore their variance is uniformly bounded 
by continuity in p and K. Thus the expression to be bounded vanishes by a law of 
large numbers. □ 

Proof of Lemma \3.1C\ The absolute value in the given equation is obviously con¬ 
tinuous in p and K. Using the fact that Aq k ^f a is closed and v(-) is continuous, it 
is straight-forward to prove that the supremum over p and K is continuous with 
respect to u. Hence for every l E N, the supremum in the hypothesis is achieved for 
some u E T d , p e Ag’t’t and 0 < K < (21 + l) d A. For these maximizers, let ip be 
a limit point as / -A oo and p p K the corresponding subsequence. Then the absolute 
value in the given formula in Lemma |3.10| is bounded from above by 


(3.13) 


v 

Fk 


■<P 


+ 


<f> 


( U > (2Z+l) d EzGAi ~ ^*( U ’ (2Z+l) d EzgA; pT”)) 


where we have used (2.4), \2.2\ and ( 3.11| ). By continuity and since p p K converges 
to <p as l -A oo, the first and last terms vanish in the limit. Now by definition of <f> 
and (2.5), we write 

<p = <f>(u,R(u,<p)) = ®(u, E m M( v{u * +qo ) ). 


Since <F(r, •) is Lipschitz continuous uniformly in u E T d , the second term in (3.13) 
is bounded from above by a multiple of 




M ( 


v(u)+q 0 


(21 + iy 


By the triangle inequality and since p E Aq ^ ,<5 a is of the form p z = q z + v z with 
\v z — v(u) | < <t, this term is bounded from above by 

L 


E M (U 


zeAi 


1 






(9) 


E„ 


M( 


v(n)+g 0 > 


1 


(2k + l) d 


E M ( 


v(n)+g^ 


zEBi 


+ Ca 
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for all small enough a > 0. Hence it is enough to consider 

1 


sup 




M( , * ) 

\ v(u)+q 0 ) 


(2k + l) d 


E M ( 


<£ 

v(u)+q z 


zE Afc 

For each q , the absolute value in this formula is bounded from above by 

n— 1 

"" L \v\u)-\-qo/ 

3=1 
n— 1 


: 1 


j =1 


+ 


n—1 


lb - _L 

, 2t TIF ^ E ("(***) ■ 

V ' 2GA fc i=0 

By the Lipschitz-continuity of M and since n < C7) _1 , this expression is bounded 
from above by C(S + aS _1 ) for all q G Aq^. Taking the limit, first for <a -A 0 and 
then for S -> 0, finishes the proof of the one block estimate. □ 

Remark 3.11. The above estimates correspond to the approximation of the inte¬ 
gral with respect to the measure m by an integral of simple functions. 

3.2. Two blocks estimate. Similarly, we can prove the two blocks estimate, which 
shows that the difference between averages over boxes of size l and eN is negligible 
in the limit. 

Lemma 3.12 (Two Blocks Estimate). It holds that 


lim sup lim sup lim sup sup E p 

l —yoo e—>-0 N—t oo |y|<eV 

m-almost surely. 


[ ik E H( x + v) ~vt N (x)\dt 

^0 rd 


We will skip the proof here and just note that it follows in a similar manner like 
the one block estimate — the main difference lying in restricting the process to two 
boxes of size l (that are at most eN apart) instead of just one. For details, the 
reader may refer to nang. 

4. Outlook and Examples: numerical simulations 

In this section, we present examples for chemotactic motion of particles in a slowly 
varying chemical attractive environment in a slightly more general setting than our 
theory has so far addressed. This is done in order to get further insight into the 
behavior of such particle systems. Consider the (smooth) function $ : T d -A (0, oo). 
Assume that the attractive chemical environment ((x) on the lattice is given 
by independent Poisson random variables with parameter at site x G T^. As 

before, let m denote the law of this random environment. One might think of this 
distribution as resulting from independent diffusion of attractive chemical molecules 
on a time-scale much faster than the typical time-scale of the diffusive motion of 
cells so that it is effectively time independent with respect to the random motion of 
the cells. This assumption on the time-scales is biologically reasonable and helpful 
also for theoretical reasons as outlined in the introduction. We suppose that the 


cells motion is described by a Markov process with generator given by (2.1), where 
the environment p N depends on the number of attractive chemical molecules £ via 

Px = /(CO)) for all x€T d N . 
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Thus / : N -A (0, oo) provides a microscopic description of the chemical bias. 
For simplicity, let us suppose for the rate g(rj) = rj , i.e. that the cells diffuse 
independently in the random environment p N . Note that the “random part” 
Px — E. m[f(C(x))] consists of independent but not stationary random variables, 
so our theory given before does not apply directly in this case. However, on local 
boxes the density is approximately constant, say 'd('u), and the environment 

consists of almost identically and independently distributed particles, so that we 
expect that with techniques as given above a hydrodynamic limit can be derived in 
this case, too. To identify this limit, consider M given as in (2.4). Since g is the 


identity, u* is a product of identical Poisson distributions with parameter <p and it 


follows M{ip) = ip. Therefore we see that R as defined in (2.5) is given by 


R(u, (/?) ®P 0 i s (i?(it)) 




L/K(*))J’ 

where the average on the right hand side is taken over a Poisson distribution with 
parameter P(u). Solving for <p yields 


P) ^Pois($(u)) 


1 


1 -1 


./(CO)). 

In other words, the effective influence of the environment is given by its harmonic 
mean. This is reminiscent of analogous formulas in (one-dimensional) stochastic 
homogenization, cf. m- Let us denote 

r 1 r 1 

/('d(R)) E Po is(^(n)) ’ 

which yields the relationship between the m icroscopic chemical bias / and its macro¬ 
scopic counterpart /. Analogous to (2.6) the equation describing the limit of the 


empirical measures is then explicitly given by 

(4.1) dtp(t,u) = A[f{d)p](t,u), (t, u) G [0, T] x T d , 


which we rewrite as 


([/(0)Vp + /'(0)pVi? 


which is of the form (1.1) with fc(p, #) = /(#) and y(p, #) = —p/'($), respectively 
%(p, $) = —/'($). From ( |4.1| ) it follows that the stationary states are given by 
multiples of 


(4.2) 


Poo (^0 — 


= E- 


Pois (i?(n)) 


/( C(x)) 


/UW) 

The respective coefficient is determined by the initial mass, which is conserved. In 
order to perform numerical simulations, let us suppose that the dependence of /(£) 
on the attractive chemical molecules is explicitly given by 

/<0 = " + IT® 

with parameters u, xo > 0 describing the relative strengths of the free diffusion and 
the chemotactic motion. Specifically, we set %o = 2 and v = 0.5. Furthermore, let 

(4.3) (tt) = 30 exp ( - 60((ui - 0.5) 2 + (u 2 - 0.5) 2 )). 


The simulations of the stochastic particle system are carried out with an algorithm 
based on random sequential updates. Here the particle densities are measured 
over balls in T 2 of radius 0.05. The results are then compared with the numerical 
solution of the corresponding limiting PDE, which is obtained by using a standard 
finite difference method. Since we only deal with the first equation of a chemotaxis 
system and because we are not in the situation where blowup in finite time is 





















16 


S. GROSSKINSKY, D. MARAHRENS, AND A. STEVENS 


expected, we do not go into further details concerning more refined algorithms for 
the PDE in this illustrating section. 

The simulations for a lattice size of N x N with N = 250 and uniform initial con¬ 
figuration Tj(x^ == 4 for all x G are shown in Figs. [l 
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Figure 1 . Cell density 
for particle system, t = 
0.0008 


Figure 2. PDE solu¬ 
tion, t = 0.0008 
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Figure 3. Cell density 
for particle system, t = 
0.004 


Figure 4. PDE solu¬ 
tion, t = 0.004 
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Figure 6 . PDE solu¬ 
tion, t = 0.01 
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for particle system, t = 
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Figure 8 . PDE solu¬ 
tion, t = 0.04 
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Figure 9. Cell density 
for particle system, t = 
0.2 


Figure 10. PDE solu¬ 
tion, t = 0.2 


The simulations for the cell density exhibit an interesting phenomenon near the 
boundary of the support of the positive density of the chemo-attractant, which it¬ 
self is not visualized in the figures. Initially a ring of lower cell density forms around 
the area where the cells finally will concentrate due to positive chemotaxis, best vis¬ 
ible as the blue ring in Figures 2 and 4. This can be explained by particles jumping 
out of the region of this ring into the domain of positive chemical attraction where 
they get stuck due to their lower motility there. If they jump out of this domain 
again, then they may not move very far away from it by random motion. This 
leads to a density profile which is a non-monotone function of the distance from the 
centre. For longer times this effect vanishes, however, and the density approaches 
the steady state, which (in the hydrodynamic limit as N -A oo) is numerically in¬ 
distinguishable from the solution of the PDE at time t = 0.2 shown in Figure [lOj 
Initially, the fluctuations of the initial condition dominate in the simulation of the 
particle system, which explains the rather large difference to the deterministic so¬ 
lution of the PDE. Averaging over several realizations of the particle system would 
lead to closer resemblance between the two. 


Finally, let us remark that it is possible to also consider models with a strong local¬ 
ization of particles. Again we assume that the chemical environment results from 
a random process. For the case considered below we do not have a proof for a 
hydrodynamic limit so far, but expect a similar outcome as before. Also, we do 
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not have stationarity in this case but convergence towards a stationary random 
variable. Now we simulate only the respective particle model, since the solution for 
the limiting PDE would require more refined numerical schemes than given above, 
because now blow-up phenomena may occur. 

Our example here is the following, let us consider a zero range process within the 
range of so-called condensation. Specifically, we consider a ZRP with a constant 
jump rate g(n) = go for all n > 0. As soon as cp > go, the partition function Z 
defined in ( |2.3| ) does not converge and the corresponding grand-canonical measure 
(2.2) does not exist. As shown in [20, 2jQ, for large enough densities the particles 
tend to concentrate on only a few sites with the lowest exit rate. However, in this 
case, this behavior is not due to the reinforced interaction of the cells and the at¬ 
tractive chemical molecules, but instead it is solely due to the stochastic properties 
of the ZRP. Such strong localization effects may relate to finite time blow-up of 
solutions of related limiting objects. 


Figs. 11 15 show the results of simulations for this model with N = 100, g{n ) = 3, 
and all other parameters given as before, i.e. uniform initial configuration rj(x) = 4, 


Xo = 2, v = 0.5 and $ given in (4.3) with final averaging over the empirical measure 
with radius 0.05. In fact, this radius can be seen explicitly at points where many 
particles have clustered in the figures below. Figure [16| shows the random density 
of the chemo-attractant in order to compare with the location of the cell aggregate. 
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Figure 16. Density of 
the chemo-attractant ft 
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If the total density in the system is high enough, we expect that eventually a fi¬ 
nite fraction of all particles will concentrate on the site with the highest number of 
chemical molecules. 


5. Concluding remarks 

In this paper we derived the first equation of a general class of chemotaxis systems 
via a hydrodynamic limit of a stochastic lattice gas. The attractive chemical envi¬ 
ronment is prescribed in our case and assumed to be random and stationary, with 
a slowly varying mean. The situation of very strong clustering of the particles on 
single lattice sites is excluded in our theory, although this phenomenon can happen 
during selforganization of chemotactic particles in certain situations, and then it 
is of important biological relevance, like for the slime mold amoebae Dictyostelium 
discoideum. It would be interesting to be able to derive the full chemotaxis system 
from a hydrodynamic limit. Technically this is more challenging and methods of 
scale separation might play a role here. 
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